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Abstract In Thailand, the forest crested lizard, Calotes 
emma consists of two subspecies, C. emma alticristatus 
Schmidt, 1925 and C. emma emma Gray, 1845. This study 
was performed to determine the genetic diversity and 
differentiation of C. emma from 16 different localities 
throughout Thailand. A total of 116 samples were 
analyzed using the mitochondrial cytochrome c oxidase 
subunit 1 (CO1). Of these, 65 and 51 of C. e. alticristatus 
and C. e. emma were classified into 23 (N1 — N23) and 21 
(S1 — S21) haplotypes, respectively. There was no shared 
haplotype between subspecies or between different 
populations within each subspecies. These haplotypes 
were classified into four (north-A to north-D) and three 
(south-A to south-C) haplogroups of C. e. alticristatus and 
C. e. emma, respectively. Phylogenetic analyses retrieved 
four lineages (classified as I to IV). Lineages I and II 
contained the four haplogroups of C. e. alticristatus, 
whereas lineage III and IV contained three haplogroups 
of C. e. emma. These two subspecies live separately in 
different climate zones, i.e. C. e. alticristatus is found in 
an equatorial winter dry climate, whereas C. e. emma 
inhabits areas with an equatorial monsoonal climate. 
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1. Introduction 


The agamid genus Calotes Cuvier, 1817 currently comprises 
23 species distributed from eastern Iran and south China 
to Sumatra, Indonesia (Vindum et al., 2003; Zug et al., 2006; 
Krishnan, 2008; Hartmann et al, 2013) with 11 species found in 
Southeast Asia (Das, 2010). Certain species have also become 
established outside of their native range (Spawls et al., 2018). In 
Thailand, C. versicolor (Daudin, 1802), C. mystaceus Buméril & 
Bibron, 1837, and C. emma Gray, 1845 occur (Chan-ard, et al., 
2015). Calotes emma has a wide distribution that ranges from 
India through Burma, southern China, Lao PDR, Cambodia, 
Vietnam and Thailand to Peninsular Malaysia (Chan-ard et al, 
2015). 

Two subspecies of C. emma, namely C. e. alticristatus Schmidt, 
1925, and C. e. emma Gray, 1845 are found in Southeast Asia (Sen 
et al., 2009; Das, 2010). Calotes e. alticristatus is distributed from 
south China (Yunnan) to Thailand, Lao PDR, Cambodia and 
Vietnam. Whereas C. e. emma occurs exclusively in Peninsular 
Thailand, Myanmar and Malaysia (Cox et al., 1998). These 
subspecies can be morphologically differentiated from each 
other by slight differences in size of the body and the spine 
above the eyes (Sen et al., 2009). The geographical divide between 
these two subspecies in Thailand was reported to occur in 
peninsular Thailand. However, the actual distribution of these 
two subspecies in Thailand and their genetic distinctiveness is 
still unclear. There is a report of C. e. alticristatus from Malaysia, 
which was identified by certain specific characters, ie. an 
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obliquely curved skin fold in front of the shoulder together 
with a small spine above the eye (Sen et al, 2009). This reveals 
that morphology alone may not be enough for subspecies 
identification or that distributional data are inaccurate. Thus 
genetic analysis should be used as an additional tool in support 
of subspecies identification. 

The patterns of genetic variation and distribution among 
populations are the results of the interactions of gene flow, drift, 
and selection processes unequally constrained by environmental 
factors across distribution ranges (Eckert et al., 2008). Different 
environments experienced by central and peripheral populations 
often results in dissimilar selection regimes thereby increasing 
genetic differentiation among populations (Conner and Hartl, 
2004). Climate is an important driver of range fragmentation 
that segregates continuously distributed populations into 
smaller, geographically isolated groups (Cheptou et al., 2017). 
Mainland Southeast Asia experiences diverse tropical and sub- 
tropical climates with highly seasonal fluctuations due to the 
monsoon and the peninsular character of the south with ocean 
on both sides. Local and regional climatic fluctuations and low 
geographic gene flow could induce severe divergent selection 
processes causing genetic differentiation between populations 
in different parts of the region. 

Genetic investigations of the genus Calotes in Southeast 
Asia are still very limited and there have been no studies 
investigating the genetic diversity within C. emma. To date, 
there is only one karyotype study comparing C. emma, C. 
versicolor and C. mystaceus in Thailand (Kritpetcharat et al., 1999). 
A population level study performed on C. versicolor from China 
and Vietnam provided evidence for high intra-population 
genetic diversity and high genetic differentiation between 
populations (Huang et al, 2013). Lastly, the genetic diversity and 
structure of C. mystaceus from Thailand and Lao PDR were 
examined based on a partial sequence of the mitochondrial CO1 
gene (Saijuntha et al, 2017). Thus, this is the first study to address 
the genetic diversity, population structure, and conservation 
measures for C. emma, across Thailand. We collected samples 
from 16 localities in northeast, central, east, west and south of 
Thailand for genetic diversity analysis, as well as for genetic 
differentiation investigations, to address the goals we mentioned 
above. 


2. Materials and Methods 


2.1. Specimen collection and DNA extraction A total of 65 
specimens of C. e. alticristatus and 51 of C. e. emma were collected 
from 16 localities in Thailand (Table 1). The lizards were caught 
using the fishing pole method (Bennett, 1999; Saijuntha et al., 
2017), then identified as C. emma based on the morphological 
characters of (i) a single spine above each eye, and (ii) twin 
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spines above the tympanum (Das, 2010). Their tails were cut 
around 5 mm from the end and soaked in 80% alcohol until 
used. After specimen collection, the lizards were photographed 
and then released back into their natural habitat. The total 
DNA was individually extracted from the tail samples using 
EZN.A. Tissue DNA kit (Omega bio-tek, USA) following the 
manufacturer's protocol. 


2.2. DNA amplification and sequencing A partial region 
of the mitochondrial cytochrome c oxidase subunit 1 (CO1) 
gene was amplified with the primers designed for C. versicolor, 
namely L5037 (5-GAG TAG ACC CAG GAA CCR AAG 
TTC-3’) and H6448 (5-GTA TAC CGG CTA ATC CAA 
GCA TGT-3’) (Huang et al, 2013). Standard polymerase chain 
reaction (PCR) was performed in 25 pL reactions, including 
approximately 100 ng of template DNA, 1 uL of each primer 
(each 10 pmol/L), 2.5 uL of 10x Ex-Taq buffer (Mg plus), 2 pL 
dNTPs (each 2.5 mmol/L), 0.125 uL of Ex-Taq DNA polymerase 
(5 U/L), and 19.375 pL of deionized water. PCR was conducted 
under the following conditions: an initial denaturing step 
at 95°C for 4 min; 35 cycles of denaturing at 94°C for 40 s, 
annealing at 65°C for 40 s, and extending at 72°C for 60 s; and 
a final extending step of 72°C for 8 min. PCR products were 
separated by electrophoresis on 1% agarose gels, visualized with 
GelRed™ Nucleic Acid Gel Stain (Biotium, Inc, Hayward, CA). 
The amplified band was cut and purified by using EZ.N.A” 
Gel Extraction kit (Omega bio-tek, USA). The purified PCR 
products were sequenced at 1" BASE Laboratories, Malaysia. All 
new sequences were deposited in GenBank with the accession 
numbers MK876009 to MK876124. 


2.3. Data analyses All sequences generated in this study 
were aligned using the ClustalW program (Larkin et al, 2007) 
and manually checked in the BioEdit program (Hall, 1999). 
Haplotype data were generated using the DnaSp v5 program 
(Librado and Rozas, 2009). A minimum spanning haplotype 
network was constructed in the Network 5.0.1.1 program based 
on the median-joining network (Bandelt et al., 1999). Molecular 
variation and neutrality test were calculated using Arlequin 
version 3.5.1.3 (Excoffier and Lischer, 2010). Genetic distance 
(p-distance) was calculated using MEGA X (Kumar et al, 2018). 
Isolation by distance (IBD) based on Mantel’s randomization test 
(Mantel, 1967) was performed in the R program (R Core Team, 
2013) ade4 package (Bougeard and Dray, 2018). We performed a 
suite of analyses to infer the phylogenetic relationship between 
C. e. alticristatus and C. e. emma from different populations 
in Thailand, as well as the four sequences of C. emma from 
Myanmar (MG935448-MG935451) available in GenBank. 
For Bayesian Inference (BI) and Mamimum Likelihood (ML) 
analysis, the best-fitting substitution model for the CO1 data 
set selected by MrModeltest ver 2.2 (Nylander, 2008) was 
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Table 1 Details of sampling localities of C. emma collected in this study. Af 


Climate zone Latitude/longitude 


Region 


Province 


District 


Sub-district 


Number 


Subspecies 


Code 
SRI 
CMI 
PLK 


14°28’22.2276” N, 101°10'1.6284” E 


Aw 
Aw 


Central 


North 


Saraburi 


Kaeng Khoi 
Fang 


Tha Maprang 


12 


C. e. alticristatus 


19°59’48.4836” N, 99°8'39.7284” E 


16°58’31.1016” N 101°2’55.8420” E 


16°44’35.8980” N, 101°1’5.3580” E 


Aw 
Aw 


15°52’20.0928” N, 101°30’44.3124” E 
17°3'13.7916” N, 101°45’1.2600” E 


Aw 
Aw 


16°29’30.2208” N, 104°22’36.9876” E 


14°30’36.8964” N, 101°55’51.2544” E 


Aw 
Aw 


North 
North 


Northeast 
Northeast 
Northeast 
Northeast 


West 


Chiang Mai 


Pong Nam Ron 


C. e. alticristatus 


Phitsanulok 


Nakhon Thai 
Khao Kho 


Ban Yaeng 


3 


C. e. alticristatus 


Phetchabun 


Khaem Son 


C. e. alticristatus 


PNB 


Chaiyaphum 


Nong Bua Rawe 


C. e. alticristatus Wang Takhe 
Puan Phu 


CPM 
LEI 


Loei 


Nong Hin 


C. e. alticristatus 


Mukdahan 


C. e. alticristatus 11 Nong Sung Nuea Nong Sung 


MDH 


NMA 


Nakhon Ratchasima 


Wang Nam Khiao 


Si Sawat 


Udom Sap 


C. e. alticristatus 


14°39'29.6064” N, 99°18'17.3160” E 


Aw 


Kanchanaburi 
Chanthaburi 


Trat 


Khao Chot 
Thap Sai 


C. e. alticristatus 


12°56’22.7508” N, 102°14’23.1252” E 


Af 


East 


Pong Nam Ron 


C. e. emma 


CTI 


12°13’47.1396” N, 102°22’53.0364” E 


9°53'40.6356” N, 98°50'1.6080” E 


Af 


Af 


Af 


East 


Laem Ngop 


Phato 


Khlong Yai 
Phra Rak 


10 


C. e. emma 


TRT 


South 
South 
South 
South 
South 


Chumphon 
Krabi 


10 


C. e. emma 


CPN 
KBI 


75530.3636” N, 99°16'5.2608” E 


Khlong Thom 


Takua Pa 


Khlong Thom Nuea 
Khuekkhak 


Nakha 
Chong 


C. e. emma 


8°37'21.2268” N, 98'14'5.0496” E 


Af 


Af 


Phangnga 


C. e. emma 


PNA 


9°27'34.5168” N, 98°30'16.1460” E 
7°33'26.9820” N, 99°46'39.4284” E 


Ranong 


Suk Samran 


3 


C. e. emma 


RNG 


Af 


Trang 


Na Yong 


C. e. emma 


TRG 
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the general time reversible with gamma distribution model 
(GTR+G). We performed a ML tree using MEGA X (Kumar 
et al, 2018) with nodal support estimated using 1000 bootstrap 
re-sampling. Whereas BI was performed using MrBayes 
version 3.1.2 (Ronquist and Huelsenbeck, 2003). The number 
of generations used in BI was 10 000 000, sampling every 100 
generations. After the first 15 000 sampled trees were discarded 
as burn-in until the average standard deviation values of the 
run dipped below 0.01 and the potential scale reduction factor 
(PSRF) value approached 1.0. The consensus tree and posterior 
probability values were calculated using the remaining 
85 000 trees. 


3. Results 


Comparison of the partial sequences of the mitochondrial CO1 
gene with a length of 1049 bp for the 116 specimens of C. emma 
confirmed that 65 and 51 specimens were genetically grouped 
according to the morphology-based subspecies classification, 
namely C. e. alticristatus and C. e. emma, respectively (Figure 1). A 
total of 134 variable sites consisting of seven singleton variable 
sites with two variants and 127 parsimony informative sites 
including 118 and nine informative sites with two and three 
variants, respectively, were observed in C. e. alticristatus. The 121 
variable sites for C. e. emma consisted of six singleton variable 
sites with two variants, and 115 parsimony informative sites 
including 106, seven, and two informative sites with two, 
three and four variants, respectively (Supplementary table 1). 
Comparison between the subspecies found that nine nucleotide 
positions were different. These variable sites can be used to 
genetically differentiate of C. e. alticristatus from C. e. emma 
(Supplementary table 1). Nucleotide and haplotype diversity 
within populations are summarized in Table 2. The number 
of haplotypes (N) within each population ranged between one 
and five haplotypes. Interestingly, no shared haplotype between 
populations was observed (Table 2). The 23 (N1-N23) and 21 
(S1-S21) haplotypes were assigned to C. e. alticristatus and C. e. 
emma respectively (Figure 1 and Supplementary Table 1). 
Genetic (p) distance ranged between 0.00290 to 0.02097 
(Table 3) and based on this together with mutational steps 
greater than 20 (Figure 1), those populations could be grouped 
as a haplogroups, thus four and three haplogroups of C. e. 
alticristatus and C. e. emma were classified, respectively (Table 
3 and Figure 1). Haplogroups of C. e. alticristatus consisted of 
north-A from Chiang Mai Province (CMI); north-B from 
Kanchanaburi (KRI); north-C from Mukdahan (MDH); 
north-D from Saraburi (SRI), Nakhon Ratchasima (NMA), 
Chaiyaphum (CPM), Phitsanulok (PLK), Phetchabun (PBN), 
and Loei (LEI). Whereas haplogroups of C. e. emma consisted of 
south-A from Ranong (RNG), Phangnga (PNA), Krabi (KBD, 
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Figure 1 Minimum spanning haplotype network (A) correlating to a map of collection sites and related to climate zones in Southeast Asia 
(B). Calotes emma alticristatus and C. emma emma were classified into N1 to N23 and S1 to S21 haplotypes, respectively. Haplotypes were 
grouped into north-A to north-D and south-A to south-C for C. e. alticristatus and C. e. emma, respectively. Climate classification is based 
on Kottek et al. (2006). Af = Equatorial fully humid, Am = Equatorial monsoonal, Aw = Equatorial winter dry, Csc = Warm temperate: 
dry and cool summer. CMI (Chiang Mai), PLK (Phitsanulok), and PBN (Phetchabun) north region; LEI (Loei), CPM (Chaiyaphum), NMA 
(Nakhon Ratchasima), and MDH (Mukdahan) northeast region; SRI (Saraburi) central region; KRI (Kanchanaburi) west region; CTI 
(Chanthaburi) and TRT (Trat) east region; CPN (Chumphon), RNG (Ranong), PNA (Phangnga), KBI (Krabi), and TRG (Trang) south 
region; MYM is where GenBank deposited sequences from Myanmar were collected. 


Chumphon (CPN); south-B from Chanthaburi (CTI) and Trat 


(TRT); south-C from Trang (TRG) (Figure 1 and Table 3). The 0.8 
north-A haplogroup of CMI was the most genetically distinct > 
with p-distance varying from 0.07820 and 0.11672. Whereas the Ë 0.6 
other species which was used as out-group, C. versicolor, had the š 
highest p-distance, the value from C. emma varied from 0.16473 z 0.4 
and 0.19750 (Table 3). 2 

Z 0.2 


An isolation-by-distance (IBD) test found that genetic 


distance and Euclidean distance between different populations 


of C. e. alticristatus and C. e. emma was significantly correlated 
with a value of R° = 0.6431 (P < 0.05) (Figure 2). The phylogenetic 


analyses of C. emma based on CO1 sequences using C. versicolor 


0 200 400 600 800 1000 1200 1400 
Euclidean distance (km) 


as out-group demonstrated four well-supported lineages (Figure Figure 2 Isolation by distance (IBD) graph plot based on Mantel’s 


3). Lineage I contained haplogroup north-A of C. e. alticristatus, randomization test between genetic distance against Euclidean 


which formed a distinct clade. Lineage II contained three distance (km) Joe 116 individual samples of Calotes emma from 16 
haplogroups, north-B, north-C and north-D of C. e. alticristatus, eae ee uae aca 

The other two lineages were of C. e emma, comprising lineage III : . 

comprising two haplogroups, south-A and south-B, and lineage 4. Discussion 
IV comprising haplogroup south-C. The four sequences of C. 
emma from Myanmar retrieved from GenBank were resolved This study is the first report of a genetic investigation of 
in a sister-group relationship with haplogroup south-B of C. e. C. emma, as well as the genetic differentiation between the 


emma from Thailand (Figure 3). subspecies, C. e. alticristatus and C. e. emma within Thailand. 
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Table 2 Summary statistics of molecular variation and neutrality analysis in 16 populations/localities of C. emma in Thailand. 


Population code N S H Uh Hd Nd 


C. emma alticristatus 


SRI 12 2 2 2 0.167 + 0.134 0.0003 + 0.0003 
CMI 3 1 2 2 0.667 + 0.314 0.0064 + 0.0003 
PLK 3 0 1 1 0.000 + 0.000 0.0000 + 0.0000 
PNB Ti 10 4 4 0.810 + 0.130 0.0045 + 0.0007 
CPM 10 2 3 3 0.511 + 0.164 0.0005 + 0.0002 
LEI 5 0 il 1 0.000 + 0.000 0.0000 + 0.0000 
MDH 11 2 3 3 0.345 + 0.172 0.0004 + 0.0002 
NMA 5 5 4 4 0.900 + 0.161 0.0027 + 0.0007 
KRI 9 3 3 3 0.556 + 0.165 0.0009 + 0.0003 
All populations 65 134 23 23 0.929 + 0.015 0.0256 + 0.0030 
C. emma emma 

CTI y 8 3 3 0.667 + 0.160 0.0043 + 0.0009 
TRT 10 0 1 1 0.000 + 0.000 0.0000 + 0.0000 
CPN 10 4 4 4 0.778 + 0.091 0.0017 + 0.0003 
KBI 4 3 3 3 0.833 + 0.222 0.0016 + 0.0006 
PNA 8 4 3 3 0.750 + 0.096 0.0019 + 0.0004 
RNG 3 1 2 2 0.667 + 0.314 0.0006 + 0.0003 
TRG 9 15 5 5 0.806 + 0120 0.0064 + 0.010 

All populations 51 121 21 21 0.939 + 0.018 0.0381 + 0.0035 


Notes: N = number of specimens, S = segregation sites, H = number of observed haplotypes, Uh = unique haplotype, Hd = haplotype diversity, 


Nd = nucleotide diversity. 


Our results revealed high genetic differences between intra- 
subspecies populations and no shared haplotypes were 
observed. The habitats of the two subspecies were separated in 
regions corresponding with different climate zones, with C. e. 
alticristatus inhabiting the equatorial winter dry climatic zone, 
and C. e. emma living in the equatorial monsoonal climatic 
zone. We found evidence for the occurrence of four lineages 
of C. emma by phylogenetic tree analyses. These lineages 
could be subdivided into seven haplogroups, which were 
classified relative to the geographical localities where they were 
found. Moreover, the genetic distances between haplogroups 
were significantly related to Euclidean distance separating 
haplogroup populations. This could suggest that there may be 
geographic barriers impacting female dispersal. 

Based on the very few discrete characters used to 
differentiate the subspecies of C. emma, many reports have 
classified these taxa at the species level only (Teynié, 2004; 
Majumder and Agarwala, 2015; Hawkeswood and Sommung, 
2017). However, there is a report of specimens morphologically 
similar to C. e. alticristatus in Jeram Linang (Lenang) Forest 
Recreation Park, Kelantan, Malaysia (Sen et al., 2009), where 


C. e. emma is known to occur. However, this could be a 


misidentification. There is also a report that C. emma does not 
occur in the west of Thailand, where C. chincollium has been 
found (Vindum 2003). This is the sister species of C. emma (Zug et 
al, 2006), and does not possess a large postorbital spine (Vindum 
et al, 2003). However, our study reveals that the specimens from 
Kanchanaburi province in the west of Thailand are genetically 
similar and cluster closely with C. e. alticristatus from the central, 
northeastern and northern areas of Thailand. Thus, the CO1 
sequence is very useful for further molecular genotyping 
and can be used to more confidentially identify/differentiate 
subspecies of C. emma throughout their distribution areas in the 
Asia-Pacific region. 

Our results provide substantial evidence that geographical 
areas relating to climate zones divide C. emma into two 
distinct lineages corresponding to the subspecies classification 
in Thailand. We found that all C. e. alticristatus existed in the 
equatorial winter dry climate zone in the north, west, and 
northeast, whereas C. e. emma was found in the equatorial 
monsoonal climate zone in east and peninsular Myanmar- 
Thailand. A population from Trang Province (TRG), the 
southernmost area in this study, was also very distinct from the 


other populations of C. e. emma and classified as another distinct 
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lineage. Interestingly a population from Chiang Mai Province 
(CMD), the northernmost area in this study, was located close to 
the warm temperate: dry and cool summer climate zone and 
was genetically very distinct from the other populations. This is 
in concordance with a previous report on C. mystaceus, which 
was also collected from the same area and was genetically very 
distinct from the other populations (Saijuntha et al, 2017). This 
finding suggests that the distribution of the two subspecies of 
C. emma is related to climate zone, as also observed in other 
organisms such as the parasitoid wasp Leptopilina boulardi 
Barbotin et al. from different climatic zones in Iran (Seyahooei 
et al., 2011). However, previous reports on the geographical 
distribution of C. emma, as well as the two subspecies, need to 
be carefully investigated by molecular analysis in order to 
determine whether the subspecies occur sympatrically. If they 
are not sympatrically distributed, then the exact distribution 
of each subspecies, as well as their genetically differentiated 
populations, needs to be explored. 
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